import java.util.Arrays;

public class FFT {
	private static final int N = 1024 * 1024;
	private static final double[] SIN = new double[N];
	private static final double[] COS = new double[N];

	static {
		if (Integer.highestOneBit(N) != N)
			throw new IllegalArgumentException("N is not a power of 2");
		for (int i = 0; i < N; i++) {
			double kth = i * -Math.PI / N;
			SIN[i] = Math.sin(kth);
			COS[i] = Math.cos(kth);
		}
	}

	private static void mult(double[] a, int k, double r, double i) {
		double ar = a[k], ai = a[k + 1];
		a[k] = ar * r - ai * i;
		a[k + 1] = ar * i + ai * r;
	}

	private static double[] fft(double[] x, int offset, int stride, int count) {
		if (count == 2)
			return new double[]{x[offset], x[offset + 1]};
		double[] q = fft(x, offset, stride * 2, count >> 1);
		double[] r = fft(x, offset + stride, stride * 2, count >> 1);
		double[] y = new double[count];
		count >>= 1;
		for (int i = 0, j = N * 2 / count, k = 0; k < count; k++, i += j) {
			mult(r, k, COS[i], SIN[i]);
			y[k] = q[k] + r[k];
			y[k + count] = q[k] - r[k];
			k++;
			y[k] = q[k] + r[k];
			y[k + count] = q[k] - r[k];
		}
		return y;
	}

	public static double[] fft(double[] ds) {
		if (Integer.highestOneBit(ds.length) != ds.length)
			throw new IllegalArgumentException("x.length is not a power of 2");
		return fft(ds, 0, 2, ds.length);
	}

	public static double[] ifft(double[] ds) {
		if (Integer.highestOneBit(ds.length) != ds.length)
			throw new IllegalArgumentException("x.length is not a power of 2");
		int n = ds.length;
		for (int i = 1; i < n; i += 2)
			ds[i] = -ds[i];
		ds = fft(ds);
		double r = 2.0 / n;
		for (int i = 0; i < n; i += 2) {
			ds[i] *= r;
			ds[i + 1] *= -r;
		}
		return ds;
	}

	public static double[] cconvolve(double[] da, double[] db) {
		int n = da.length;
		if (n != db.length)
			throw new IllegalArgumentException("x.length != y.length");
		double[] a = fft(da);
		double[] b = fft(db);
		for (int i = 0; i < n; i += 2)
			mult(a, i, b[i], b[i + 1]);
		return ifft(a);
	}

	public static double[] convolve(double[] da, double[] db) {
		return cconvolve(Arrays.copyOf(da, da.length * 2), Arrays.copyOf(db, db.length * 2));
	}

	public static double warmup() {
		double[] x = new double[N * 2];
		for (int i = 0; i < N * 2; i += 2)
			x[i] = (double)i / -N + 1;
		double[] y = fft(x);
		double[] z = ifft(y);
		double[] c = cconvolve(x, x);
		double[] d = convolve(x, x);

		double r = 0;
		for (int i = 0; i < N * 2; i += 2)
			r += y[i] + z[i] + c[i] + d[i];
		return r;
	}

	public static void main(String[] args) {
		for (int i = 0; i < 5; i++)
			warmup();
		for (int j = 0; j < 5; j++) {
			long begintime = System.currentTimeMillis();
			double r = warmup();
			System.out.println(System.currentTimeMillis() - begintime);
			System.out.println(r); // 输出结果,避免计算过程被编译器优化掉
		}
	}
}
